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ABSTRACT 

o 
o i 

We present the results of Monte Carlo simulations for the dynamical evolution 
of star clusters containing two stellar populations with individual masses mi and 
m 2 > m i, and total masses M\ and Mi < Mi. We use both King and Plummer model 
initial conditions and we perform simulations for a wide range of individual and total 
fT} ■ mass ratios, mil "mi and Mi/ Mi. We ignore the effects of binaries, stellar evolution, 

and the galactic tidal field. The simulations use N = 10 5 stars and follow the evolution 
fN| . of the clusters until core collapse. We find that the departure from energy equipartition 

, in the core follows approximately the theoretical predictions of Spitzer (1969) and 

I/" - ) | Lightman & Fall (1978), and we suggest a more exact condition that is based on our 

results. We find good agreement with previous results obtained by other methods 
regarding several important features of the evolution, including the pre-collapse 
distribution of heavier stars, the time scale on which equipartition is approached, and 
the extent to which core collapse is accelerated by a small subpopulation of heavier 
stars. We briefly discuss the possible implications of our results for the dynamical 
q . evolution of primordial black holes and neutron stars in globular clusters. 

CO 

Subject headings: clusters: globular — celestial mechanics, stellar dynamics - Monte 
Carlo: dynamical evolution 



1. Introduction 

Remarkable advances have been made over the last three decades in our understanding of 
globular cluster dynamics (see, e.g., Meylan & Heggie 1997 for a recent review). The simple 
case of a two-component cluster is traditionally regarded as the second level of sophistication 
and therefore a logical challenge for new methods that have tackled the single-component case. 



Present address: Laboratoire de Physique de la Matiere Condensee, Ecole Polytechnique, 91128 Palaiseau, France; 
email: ww@pmc.polyteclinique.fr. 

2 6-218M MIT, 77 Massachusetts Ave, Cambridge, MA 02139; email: kjoshi@mit.edu. 

3 6-201 MIT, 77 Massachusetts Ave, Cambridge, MA 02139; email: rasio@mit.edu. 

4 Alfred P. Sloan Research Fellow. 



- 2 - 



Two-component clusters were originally examined because they better resemble real clusters, 
which contain a continuous spectrum of masses. While somewhat more realistic in this regard, 
clusters with only two mass components still represent a simplification with respect to real 
clusters. It has been suggested recently, however, that for a range of configurations in mass types 
and the relative size of the two populations, two-component clusters can resemble real clusters 
that are mostly comprised of compact objects and main-sequence stars (Kim, Lee, &; Goodman 
1998). Similarly, clusters containing both single main-sequence stars and primordial binaries can 
be modeled, in first approximation, as two-component systems, although dynamical interactions 
involving binaries are expected to play an important role for these systems (Gao et al. 1991). 
Perhaps the best reason to examine a simplified model of any stellar system, however, is to obtain 
a more profound understanding of individual physical processes. 

Much of the discourse regarding two-component systems has focused upon the following 
questions. First, for what configurations of the cluster is dynamical equilibrium precluded (i.e., 
the system is not stable on dynamical time scales)? Second, for what configurations of the 
cluster is thermal equilibrium precluded (i.e., equipartition of kinetic energies between each 
component is not allowed)? Both questions originate from an analysis by Spitzer (1969), in which 
he noticed that simultaneous thermal and dynamical equilibrium is impossible for some clusters. 
In particular, the heavier stars sink into the center as they lose kinetic energy to the lighter stars 
during the approach to equipartition. If equipartition is not attained, then the heavier stars will 
continue sinking until their self-gravity dominates the cluster's potential in the core. Shortly 
thereafter, the heavier component will undergo a gravothermal collapse, forming a small dense core 
comprised mainly of the heavier stars ( |Spitzer 1969 ). Refinements to this analysis have obtained 



similar constraints upon the configurations of two-component clusters in dynamical and thermal 
equilibrium (Lightman & Fall 1978). 



Several methods have been used to address questions about dynamical and thermal equilibrium 
in two-component systems. These include the construction and study of one-parameter families of 
models in dynamical equilibrium flKondrat'ev &; Ozernoy 1982 ; Katz & Taff 1983), Monte Carlo 



approaches to the numerical integration of the Fokker-Planck equation ( |Spitzer Hart 197l| ), 
direct integration of the Fokker-Planck equation in phase space ( llnagaki fc Wiyanto 1981 |K lmJ 



Lee, fc Goodman 1998 ), and also direct iV-body simulations (Portegies Zwart & McMillan 2000). 
The majority of work using any one of these methods has been undertaken at least partly in order 
to confirm Spitzer's conclusion flYoshizawa et al. 1978 ) or refute it (Merritt 1981). 



Dynamical equilibrium is attained and maintained on time scales that are very short compared 
to the amount of time needed for relaxation or equipartition. A so-called "equilibrium model" 
(i.e., whose phase-space distribution function satisfies the equation for hydrostatic equilibrium) 
therefore resembles a possible stage or snapshot in the evolution of a dynamically stable cluster. 
It is interesting to construct a parametrized family of equilibrium models, for which equipartition 
is either assumed or the temperature ratio allowed to vary, in order to determine under what 
conditions the dynamical equilibrium becomes impossible. In the majority of previous work, the 
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distribution functions of such families take the form of lowered Maxwellians or spatially-truncated 
isothermal spheres. 

Yoshizawa et al. (1978) examined a linear series of equilibrium models of two-component 
isothermal spheres with reflecting walls and found that turning points of the total energy (stability 
limits) were positive for some clusters. In such cases, a cluster was not self-bounded and would 
dissociate if the walls were removed. These configurations, which were interpreted to preclude 
dynamical equilibrium, occurred under conditions corresponding closely to the ones Spitzer 
proposed. Katz & Taff (1983) examined the limits of stability for a linear series of equilibrium 
models of self-bounded two-component clusters with a lowered-Maxwellian velocity distribution. 
They found that the number of possible configurations for clusters in dynamical and thermal 
equilibrium were diminished dramatically in the regime that Spitzer proposed. Kondrat'ev & 
Ozernoy (1982) constructed a family of equilibrium models, based upon a generalization of the 
single-component King models, for which equipartition of kinetic energies was generally impossible. 
By contrast, Merritt (1981) described a one-parameter family of equilibrium models for which 
equipartition in the core is possible for all cluster configurations within Spitzer's unstable regime. 
It seems clear, however, that these models violate an important assumption of Spitzer's analysis, 
and that they are highly unrealistic ( [Merritt 1981 ). 

Few studies that examine the evolution of two-component systems have been undertaken for 
a wide range of cluster configurations and using the other methods mentioned above. Among the 
more notable efforts are the Monte Carlo calculations of Spitzer & Hart (1971; later extended to 
three-component systems by Spitzer & Shull 1975), which were interpreted to partially confirm 
Spitzer's original analysis, and the direct Fokker-Planck integrations of Inagaki & Wiyanto (1984; 
see also Inagaki 1985), which also partially support Spitzer's claim, although by examining a 
limited range of clusters that do not totally satisfy his assumptions. The former study finds that 
a single model, which belongs to the unstable regime in which equipartition is precluded, develops 
a collapsing subsystem of heavier stars. Both studies find that central equipartition of kinetic 
energies is never attained throughout the evolution of some clusters. 

In this paper we present the results of calculations used to model the evolution of two- 
component clusters with a wide range of configurations. We examine several features of the 
evolution, including the tendency toward equipartition, the evolution of mass densities in the core, 
and the tendency for core collapse to be accelerated by the presence of a small subpopulation of 
heavier stars. Our aim is partly to assess the accuracy of Spitzer's analysis and also the refinement 
supplied by Lightman & Fall (1978). We have used a new Monte Carlo code for modeling the 
evolution as a sequence of equilibrium models whose velocities are perturbed according to the 
average effect of long-range stellar encounters (Joshi, Rasio, & Portegies Zwart 2000a). Our initial 
two-component systems are isolated King and Plummer models with a subpopulation of heavier 
stars. 

Our paper is organized as follows. In §^ we review the theoretical arguments that suggest 
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conditions under which simultaneous dynamical and thermal equilibrium are not allowed. In §|3] 
we discuss our numerical method and the quantities we calculate, and in §|3] we present our main 
results. In the last section, §[5], we discuss some additional comparisons between our findings and 
those of other studies, as well as some astrophysical implications of our results for subpopulations 
of compact objects in globular clusters. 



2. Equipartition of Kinetic Energies in Two- Component Clusters 

According to Spitzer (1969), simultaneous dynamical and thermal equilibrium is impossible 
for two-component star clusters within a certain range of configurations. Let us consider a cluster 
with stars of two masses, mi and m 2 -, where m 2 > mi. Moreover, let M 2 and M\ be the total 
mass in each component. Spitzer assumed that M 2 <C Mi, which is normally the case for real 
clusters. He concluded that for certain values of the ratios M2/M1 and m 2 /mi, equipartition will 
not be attained as the heavier and lighter stars exchange kinetic energy, and hence the heavier 
stars will sink very far into the center. Moreover, the heat exchange with lighter stars promotes 
them to higher orbits, so that eventually insufficient numbers will remain in the center to conduct 
heat rapidly away from the heavier stars. If the mass-stratification proceeds far enough, then the 
self-gravity of the heavier stars will dominate the potential in the core, and this subsystem will 
undergo gravothermal collapse. The result is a very dense core comprised exclusively of heavier 
stars. 

In his analysis, Spitzer begins by assuming global equipartition. Global equipartition is 
not realistic, however, because the relaxation and equipartition times vary greatly throughout 
the cluster, becoming longer than the age of the universe in the outer halo. In fact, we expect 
equipartition only in the inner region where relaxation times are shortest. His discussion is mostly 
unchanged by this, so long as we confine its relevance to processes in this inner region. We shall 
reproduce here only the main strategy of his argument and its conclusions. Let v^i and t>^ 2 
represent the mean-square velocities of stars in each component. As mentioned, The assumption 
of equipartition implies that the temperature ratio £ is equal to unity, 



(l/2)m 2 < 2 
(l/2)mn4i 



Spitzer then applies a component-wise virial theorem, which for the heavier component gives 

2 _ (2/5)GM 2 G r P2M l (r)dV 

v m2 — I r -rz- I - , {*) 



r 2 M 2 Jo 

where the first term represents the gravitational self-binding energy of the heavier stars, r 2 
is the virial radius of the heavier stars, p 2 is the density of the heavier stars at a distance r 
from the center, Mi(r) is the total mass of the heavier stars contained within the radius r, 
and dV is the volume element. Spitzer assumes that, since M 2 -C Mi, the second term in the 
corresponding equation for be ignored, and also that, since m 2 3> mi, the heavier stars 
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will be concentrated in the center of the system and Mi(r) may be approximated by 4-7r / oi(0)r 3 /3 
in equation (|2|). Spitzer & Hart (1971) noticed that the second assumption often does not hold 
strictly. In particular, even far into the evolution many heavier stars can still reside well outside 
of the core. Merritt (1981) constructed equilibrium models that violate this assumption by 
great amounts and discovered that equipartition is possible for some (admittedly unrealistic) 
configurations which can be realized for all values of m^jrax and M 2 jM\. 

Under Spitzer's assumptions, after a series of manipulations one obtains the following 
expression for the quantity S, as a direct consequence of equation ([[]), 

/M 2 \ /m 2 \ 3 / 2 = ( Phl / Ph2 ) 1/2 
*-\Mj\mJ (l + W/%2) 3 / 2 ' K1 

where phi and ph2 are the densities within the half- mass radii for each component, and where a 
depends on the distribution of mass within the cluster. In particular, if we denote by r m 2 and r^2 
the rms and half-mass radius for the heavier component, respectively, then a is given by 

a= 5M0)/r^y (4) 
4p/ri \r h2 J 

Spitzer estimates a value of 5.6 for a. Merritt (1981) contests the value 3.5 assigned by Spitzer to 
pi(0)/phi for polytropes with n between 3 and 5, but finds that it corresponds approximately to its 
minimum value, which, as we shall see, does not change Spitzer's conclusion. In particular, since 
the RHS of equation (||) has a maximum value 0.38a -1 / 2 with respect to variation in Phi/ph2, it 
follows that equipartition is possible only if S does not exceed a critical value (3 = 0.38a -1 / 2 , 

S < (3, where = 0.16 for a = 5.6. (5) 

Spitzer suggests that for smaller values of the individual mass ratio 7712/771-1, the inequality (||) 
remains valid, except that (3 — > 1 as 7772/7711 — > 1. 

While it is useful to assess under what conditions simultaneous dynamical and thermal 
equilibrium are not expected, it is also interesting for our purposes to estimate the extent of 
departure from thermal equilibrium where dynamical equilibrium holds. To that end, let us 
maintain the condition stated in (^), while not insisting that £ = 1 in equation (|l|). By following 
the value £ through Spitzer's analysis, we find that the temperature ratio has a lower bound for a 
given value of 5, 

?>y . (o) 

Inagaki & Wiyanto (1984) performed an analysis similar to Spitzer's, except by casting 
equations (|l|) and (|2|) in terms of core values for the component-wise total mass, mean-square 
velocity, and density. They obtained a minimum difference between the core temperatures of each 
component. Letting M c2 and M c \ denote the total mass contained within the core for the heavier 
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and lighter components, respectively, and letting r c i represent the core-radius of the lighter stars, 
this difference is given by: 



1 



-m 2 (u 



1 



m2Jmm 



-m x v m2 



1 (M c2 \ 

2 \M cl ) 



2/3 



0.2mi 



GMd 
r c i 



0) 



For now, we note that the minimum ratio (given by the inequality [6]), and difference (given by 
eq. [7]) of component-wise temperatures increases with increasing values of the total mass ratio, 
and that this increase is steeper for increasing m 2 /mi (in eq. [6], recall that [5 — > 1 as m 2 /m\ — > 1). 

Lightman & Fall (1978) developed an approximate theory for the core collapse of two- 
component clusters which resembles that of Spitzer. They examined two constant-density 
isothermal spheres representing the cores of the heavier and lighter components, where the radius 
of the former is smallest. By applying a component- wise virial theorem and several simplifying 
assumptions, they find a set of four ordinary differential equations for the virial radii and total 
masses in each component. They obtain the following condition for equipartition of energies 
between the two components, where we let m = m 2 jm\ and M = M 2 /M\: 



r(m, M) = m 3 M 2 



3MW 5 ~ 

1 H 1 + -M 

2m I V 2 



-3 



< 1 



(8) 



Solutions to the differential equations exhibit a minimum temperature ratio f; m in> which they 
suggest bears the following approximate relation to T, 

Cmin — F I . (9) 

In this case also, the minimum temperature ratio increases with increasing values of the individual 
and total mass ratios. 

The majority of attempts to evaluate the theoretical predictions of Spitzer (1969) and 
Lightman &: Fall (1978) have met with moderate success. For example, Yoshizawa et al. (1978) 
obtained 0.25 for the value of f3 in the case of spatially-truncated two-component isothermal 
spheres. Nevertheless, few investigations have examined models with a comprehensive range of 
individual and total mass ratios (outside of studies based upon turning points along a sequence of 
equilibrium models). In the next two sections we discuss the methods we used and results that we 
obtained for a relatively broad survey of the parameter space determined by M 2 jM\ and m 2 jm\ . 



3. Numerical Methods and Definitions 

We used a Monte Carlo method for modeling the dynamical evolution of clusters as a sequence 
of equilibrium models subject to regular velocity perturbations. The velocity perturbations 
represent the average effect of long-range stellar encounters ( [Henon 1971 ). Our Monte Carlo code 



has been described in detail by Joshi, Rasio, & Portegies Zwart (2000a) and Joshi, Nave, & Rasio 
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(2000b). The code allows us to perform dynamical simulations for realistic clusters containing up 
to N ~ 10 5 - 10 6 stars on a parallel supercomputer. In the present study, we ignore the effects of 
binaries, stellar evolution, and the galactic tidal field. Our units are defined in Joshi et al. (2000a, 
§2.8). The unit of length is close to the virial radius of the cluster, the mass is measured in units 
of the total initial cluster mass, and the unit of time is of order the initial half-mass relaxation 
time t r h. In this paper, when reporting times in units of i r h, we have adopted the same definition 
used in previous studies of two-component clusters since Spitzer & Hart (1971), 

tA ~ GV2 m l ogl0 (0.47V)' (10) 

where M = Mi + M<i is the total cluster mass, is the initial cluster half-mass radius, and 
m = M/N is the average stellar mass. 

We undertook two sets of calculations, hereafter called sets A and B. All calculations are 
performed for a cluster containing N = 10 5 single stars. The initial model used in each calculation 
of set A was a two-component King model ( [King 196G| ). In particular, the velocities and positions 
for all stars with a mass mi were chosen according to the King model distribution function with 
dimensionless central potential Wo = 6. Although the initial King model is truncated at its finite 
tidal radius, we do not enforce that tidal boundary during the evolution, allowing the cluster 
to expand indefinitely. In each case, some fraction of the stars were then changed to a mass 
m,2 > mi according to a chosen value of the total mass ratio M%/Mi. The initial ratio of mean 
temperatures in the heavier and lighter components was therefore m<ijmi. The initial models for 
calculations in the set B were generated in a similar way, except using the Plummer distribution 
function (polytrope with n = 5; see, e.g., Binney & Tremaine 1989). Both sets of calculations 
are listed in Table |]. The set A includes calculations undertaken for a range of total mass ratios 
(M%jMi = 3 x 10~ 3 to 0.6) and a range of individual mass ratios (m^/mi = 1.5 to 6). The set B 
comprises only 9 systems, including 4 studied by Inagaki & Wiyanto (1984), all with m^jmi = 2. 
Every calculation is terminated at core collapse, measured at the instant that a number density 
of 10 8 (in our units) is attained in the core. Our results are not sensitive to the exact value of 
the core density used to terminate the calculation. However, the value of the core-collapse time 
t cc determined numerically can have a large statistical uncertainty, particularly when the core is 
dominated by a small number of heavier stars near the end of the evolution (in those cases we 
estimate that the statistical uncertainty on the values of i C cArh reported in Table 1 can be as large 
as ~ 5%). The majority of our calculations require between 10 — 20 CPU hours on an SGI/Cray 
Origin2000 supercomputer to reach core collapse. 

The range of values we consider for mil mi and M-z/Mi includes a number of astrophysically 
relevant cases. In particular, a subpopulation of neutron stars in a dense globular cluster might 
have m%lm\ ~ 2 (e.g., corresponding to mi = 1.4 M Q for an average background stellar mass 
mi = 0.7 Mq) and M2/M1 ~ 10 -3 — 10 -2 depending on the neutron star retention fraction. A 
subpopulation of black holes might have m^jvai ~ 5 — 10 and M2/M1 ~ 10~ 3 — 10~ 2 . Massive 
blue stragglers or primordial binaries containing main-sequence stars near the turn-off mass, would 
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have m 2 /m 1 ~ 2 - 3 and M 2 /M 1 ~ 1(T 3 - 10 _1 . 

For each calculation we record several quantities at each program time step (the time steps 
are proportional to a fraction of the core relaxation time; see Joshi et al. 2000a). These continuous 
measurements include the total-cluster core radius and several component-wise Lagrange radii. 
Within each of these radii we also count the number of stars and calculate the mean temperature 
and mean mass density for each component. Of particular interest are the quantities calculated 
inside the core radius, where relaxation times are shortest and where thermal equilibrium is the 
most likely to occur. We use the customary definition for the total-cluster core radius r c flSpitzei 



19871) , 

/ , >,\ 1/2 

_ / 3v m (0) 2 \ 



(11) 



\4irGp(0) J ' 

where v m (0) 2 is the mean-square velocity and p(0) the mean density of stars at the center. 
We calculate the ratio of core densities in each component (p C 2/Pcl), and also the ratio of core 
temperatures. We also calculate the minimum core temperature ratio £ m ; n that is reached after 
approximately 90% of the pre-collapse evolution in all cases. Specifically, £ m i n is calculated as 
the average temperature ratio from 90% to 95% of the core collapse time. The core collapse 
time t cc and the time t p at which the core mass densities of each component become equal (i.e., 
Pcil Pc\ = 1) are also measured. We report our main results in Table [j], and we discuss these in the 
following section. 



4. Results 

The evolution of the core temperature ratio £ is shown in Figure [I] for three calculations in the 
set A (two-component King models): namely, for S = 0.05 (top), S = 0.5 (middle), and S = 1.24 
(bottom). Figure 2 shows the core temperatures of the lighter stars (top) and the heavier stars 
(bottom) for the case S = 1.24. Several features that we expect and that have been mentioned 
already in §|2] are easily recognizable. The temperature ratio begins with a value mj/mi and 
then decreases gradually as equipartition is approached. It is clear that £ reaches a minimum 
value that is greater than 1 for the case S = 1.24, so that equipartition is clearly never attained. 
Equipartition is nearly attained for S = 0.5, and £ m j n = 1 to within 5% for S = 0.05. It is clear 
from Figure [2] that the heavier component cools initially, then maintains a constant mean kinetic 
energy, and then begins to heat prior to core collapse. At the same time, the lighter component 
steadily becomes hotter as it receives energy from the heavier component. The temperature 
ratio in the last time steps becomes very noisy because the temperatures are computed using the 
relatively few stars that remain in the core. 

The temperature ratio £ reaches a minimum value at different times with respect to core 
collapse for each of the models shown in Figure In cases where the minimum value is greater 
than 1, it sometimes appears that the gravothermal catastrophe has beaten the approach to 
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equipartition. In such cases, one may ask whether an initial model with a less concentrated spatial 
distribution and a different initial relaxation time would yield a different minimum temperature 
ratio. In fact, we find that £ m i n is robust with respect to changes in the initial value of the 
dimensionless potential Wo- The evolution of the temperature ratio for three calculations, where 
Wo = 1, 5, and 10 for S = 1 and mijm\ = 5, are shown in Figure [|. In all three cases the 
minimum temperature ratio is approximately 1.55. 

The evolution of the core mass density ratio p C 2/Pci is shown in Figure || for the same three 
clusters and in the same order. The core mass densities of each component are shown in Figure |H| 
for the model with S = 1.24. One can see clearly that, as S increases, the core mass densities 
in each component become equal sooner with respect to core collapse. That is, for larger S the 
self-gravity of the heavier stars dominates the potential in the core for a longer period prior to the 
onset of the gravothermal catastrophe. Moreover, we can see that, although the core density is 
initially dominated by the lighter stars, the heavier stars overtake and exceed the density of lighter 
stars by more than an order of magnitude prior to core collapse for S = 0.5 and S = 1.24. 

Approximate values of the minimum core temperature ratio ^ m j n are plotted using three 
symbols in the parameter space determined by MijM.\ and mijmx in Figure ^ for 37 calculations 
in the set A. Also drawn are the Spitzer and Lightman-Fall "stability boundaries," above which 
simultaneous dynamical and thermal equilibrium are supposedly prohibited (S = 0.16 and T = 1, 
respectively; cf. eqs. [5] k, [8]). Our simulations allow us to determine £ m i n with a numerical 
accuracy of about 5%. Specifically, £ m i n is calculated as the average core temperature ratio 
from 90% to 95% of the pre-collapse evolution, and this average has a standard deviation of 
approximately 0.05 in our calculations for N = 10 5 stars. Therefore, those calculations marked 
with a "□" in Figure |6| have been determined to reach equipartition within our numerical accuracy. 
One can see that the Spitzer boundary S = 0.16 is approximately respected for 777,2/7711 > 2. By 
comparison, the Lightman-Fall boundary falls well inside the range of clusters which have clearly 
not attained equipartition. In spite of this, the Lightman-Fall boundary better reproduces the 
shape of boundaries between regions of constant £ m in- A more properly-drawn Spitzer boundary 
has a similar shape, recalling that (3 — ► 1 as m^/mx — ► 1. Based upon the results shown in 
Figure |6|, we propose our own condition for equipartition, 

L=mm"z°*. (12) 



Mi J \mi 



The boundary determined by equation (|T2|) is strictly valid for 1.75 < 777.2/?™ 1 < 7, and is also 



drawn in Figure |(| For mijfrix < 1.75, equipartition is achieved for all clusters considered. For 



railra\ > 1.75, extrapolation of equation (12) seems reasonable. 

The dependence of £ m j n on S for several values of mj/mi is shown in Figure 0. Also drawn is 
the Spitzer stability boundary. These curves are broadly consistent with trends anticipated by the 
inequality (6). In particular, £ m j n increases with S, and the initial slope of this increase becomes 
larger with increasing mi/mi (again recalling that (3 — > 1 as mijvnx — * 1)- The dependence of £, 



mm 



on r for several values of 7712/771-1 is shown in Figure |8[ Also drawn is the value of £ m i n anticipated 
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by equation (Q), and the Lightman-Fall stability boundary. One can see that, while the numerical 
results are displaced from the predicted curve, they have the same curvature and display a 
similar trend. The minimum temperature difference for calculations in the set B (two-component 
Plummer models) are shown in Figure [], where they are compared with the results of Inagaki 
& Wiyanto (1984). Here, agreement is excellent except in the limit of small M2/M1, where the 
temperature of the heavier component is calculated using very few stars, so that the difference is 
characterized by a large amount of noise. As a last comparison, we note that, using a Monte Carlo 
scheme different from ours, Spitzer & Hart (1971) found that £ m ; n = 1.34 for a Plummer model 
with S = 1.24 and m2/m\ = 5, whereas for the same system we obtain £ m j n = 1.60. 

Recall that t p is the time at which the mass densities of each component become equal in 
the core. Approximate values of the time t p as a fraction of the core collapse time t cc are plotted 
using three symbols in the parameter space determined by M2/M1 and mz/mi in Figure [l^, for 37 
calculations in the set A. This plot confirms the trend anticipated by the previous examination of 
three individual clusters (Fig. 4). Namely, the amount of time — as a fraction of the core-collapse 
time — during which the heavier stars dominate the potential in the core increases with S. 
The trends appear not to respect any of the previous stability boundaries very well, but our 
condition (12) fares best. They may nevertheless shed light on the related question of whether 
a dense subsystem of heavy stars collapses independently, since the self-gravity of the heavier 
stars must dominate the potential in the core in order for this to happen. Where equipartition is 
attained, the mass segregation is retarded or stopped, so that equal mass densities may not occur 
until core collapse (i.e., t p ~ t cc ). 

All of the calculations were terminated at core-collapse, at which time the radius containing 
1% of the mass in the heavier component diminishes sharply. The time t cc is measured at the 
instant when a number density of 10 8 in our units is attained within the core (see §3). The 



variation of the core collapse time t cc with S for several values of m2/m-i is shown in Figure 11 



The trends confirm that the onset of core collapse is accelerated by the presence of a small and 



heavier subpopulation, in agreement with the findings of others (inagaki & Wiyanto 1984; Inagaki 
1985; iQuinlan 19961) . 



5. Discussion 

In this section we discuss several features of the evolution in more detail, and we mention 
some possible astrophysical applications of our results to the evolution of compact stellar remnants 
in globular clusters. 

The temperature ratio £ in Figure Q initially has the value 7112 /mi. While this is an artifact 
of the way our initial models were constructed, m^jvrtx happens also to be the most realistic value 
of £ for equilibrium models with a relatively shallow potential. In families of equilibrium models it 
is typical to find that f -> m 2 /m 1 as W -> ( [Kondrat'ev fc Ozernoy 1982] , [Katz fc Taff 1983|) . In 



- 11 - 



trials for which initial models were modified so that £ had some value other than mijm\^ a brief 
period of rapid relaxation ensued which restored the value ni2/mi. This effect has been observed 
in simpler models of the evolution calculated using other methods (Lightman & Fall 1978). 



In the core, the initial behavior of the temperature ratio is mostly determined by the 
temperature of the heavier component, while the mean kinetic energy of the lighter stars, which 
are more abundant at first, increases gradually (Fig. |2|). Spitzer (1969) suggested that the 
approach to equipartition is characterized by the exponential decay of kinetic energy in the heavier 
component, with a time constant equal to twice the so-called equipartition time, 



t eq = t rl — i 1 + , (13) 

16 m 2 V v m iJ 

where t r \ is a relaxation time for the stars of mass m\. In the case where the mean-square 
velocities of each component are initially equal, the initial equipartition time is i cq ~ t r \{m\jm2). 
(It should be noted that t eq decreases as equipartition is approached.) This characterization of 
the decline in kinetic energy of heavier stars agrees very well with our results for stars contained 
within the half- mass radius, where we assume t r \ ~ i r h, the initial half-mass relaxation time for 
the cluster as a whole. In particular, the kinetic energy of the heavier component diminishes to 
a fraction 1/e of its initial value (after subtracting its minimum value for the entire evolution) 
within 0.39t r h for 5 = 1, mijrax = 5, and within 1.3i r h for S = 0.5, mijm\ = 1.5, in good 
agreement with the theory (which predicts a time 2t cq ~ 2{m\ / m2)t r f l ) . However, we find that 
equipartition is approached on a similar time scale in the core, where the theory predicts that t eq 
should be shorter by approximately 1/5, and hence agreement is poor (the ratio of initial core and 
half-mass relaxation times for King models with Wq = 6 is approximately 1/5; see Quinlan 1996] ). 



A leveling in the temperature ratio at a minimum value greater than 1 has been observed in 
calculations using other methods as well ( |Inagaki fc Wiyanto 1984 , Lightman fc Fall 1978; ). We 



find that this leveling is approximately coincident with the heavier stars reaching their maximum 
numbers within the core. Inagaki & Wiyanto (1984) found that an increase in the core temperature 
ratio prior to or during collapse is coincident with t p , the time at which equal mass densities are 
attained in the core. While we are not able to resolve adequately the late-collapse behavior of £ 
(because our calculation loses accuracy in this regime), it is clear from Figure |2| that the heavier 
component begins to heat prior to collapse. Indeed, we find that the temperature of the heavier 
component does not begin to rise until t > t p . 

Katz & Taff (1983) examined the turning points in a linear series of equilibrium models. In 
particular, they studied a one-parameter family of self-bounded isolated equilibrium models with 
a lowered-Maxwellian velocity distribution. Turning points representing the limits of stability 
for dynamical equilibrium were obtained for several values of mijm\ and M2/M1 in terms of 
maximum possible values of the dimensionless potential k. Katz & Taff also calculated the core 
temperature ratio £ that corresponds to each maximum value of k. Since for their models £ was 
found to approach 1 for large values of k, the core temperature ratios calculated for each turning 
point represent the minimum allowed value of £ for given values of 1712/ mi and M 2 /M\. These 
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are plotted with respect to S for several values of 7772/7771 in Figure and should be compared 
with our results, shown in Figure [7|. While these numbers exhibit a similar trend with S for given 
777,2/777-1, the trend in 7772/7)71 evidently disagrees with our results, and therefore also the prediction 
of the inequality (g). It is important to bear in mind, however, that these ratios are obtained 
just for the members of this particular family, which do not necessarily represent states that can 
be obtained by real dynamical processes. Our results contradict those of Katz & Taff (1983) 
insofar as we find clusters with smaller values of £ m j n for identical values of 7772/7771 and M2/M1, 
and which appear to be stable on dynamical time scales. Katz & Taff found that the number of 
possible configurations for their models in dynamical equilibrium diminish sharply for S > 0.16. It 
is interesting to note also that the values they obtain for £ m j n appear not to diminish below 1.10 
for small values of M%/Mi. 

We concur with the findings of Spitzer & Hart (1971) that many heavier stars remain outside 
the core throughout the evolution. This is clear from plots of Lagrange radii as a function of time 
for the heavier component (Fig. 13). This casts doubt on the assumption, committed in Spitzer 's 



original analysis, that all of the heavier stars quickly become concentrated in the core (see §y). 
In particular, we find that for S = 1.24, M2/M1 = 0.111, and 7772/7771 = 5, the radius containing 
75% of the mass in the heavier component diminishes to only 50% of its initial value (and hence 
remains larger than the core radius) throughout the evolution. Nevertheless, by the onset of core 
collapse, we frequently observe for calculations with large mijra\ that no lighter stars remain in 
the core. 

Our results may have important implications for the dynamical evolution of various 
subpopulations of interesting objects in globular clusters. In particular, a subcomponent of 
primordial black holes with 7772/7771 ~ 10 is expected to remain far from energy equipartition with 
the rest of the cluster, and to evolve very quickly to core collapse on its own relaxation time 
scale. For a typical cluster IMF (initial mass function), and assuming that all black holes formed 
initially by the stellar population are retained in the cluster, we expect M 2 /M 1 ~ 10~ 3 - 10~ 2 , 
well above our stability boundary in Figure 6. Recent iV-body simulations for clusters containing 
primordial black holes indicate that, after reaching core collapse, the dense subcluster of black 
holes evaporates quickly in the background cluster potential. Three-body processes occurring in 
the post-collapse phase produce a significant number of tight black hole binaries that will coalesce 
in a few billion years, making these binaries important sources of gravitational waves for current 
ground-based laser interferometers (Portegies Zwart & McMillan 2000; see also Kulkarni, Hut, & 
McMillan 1993 and Sigurdsson & Hernquist 1993). For neutron stars, with 7772/7771 ~ 2, our results 
suggest that equipartition can be reached if the fraction of the total cluster mass in neutron stars 
is ^ 5%. This fraction is higher than would be predicted for a standard cluster IMF and neutron 
star progenitor masses (even if, in contrast to what is suggested by many observations, neutron 
stars were born without the kicks that might eject a large fraction from the cluster). However, 
many multi-mass King models and dynamically evolving Fokker-Planck models of globular clusters 
based on fits to both photometric and kinematic data suggest that much higher fractions of 



-13- 



neutron stars may be present in many clusters. For example, the recent Fokker-Planck models 
of Murphy et al. (1998) for 47 Tuc suggest that 4.6% of the total cluster mass is in the form of 
dark stellar remnants of mass 1.4M . With such a high mass fraction, it is possible that the 
neutron stars in 47 Tuc may remain out of energy equipartition with the rest of the cluster. More 
sophisticated dynamical simulations taking into account the full mass spectrum of the cluster, 
stellar evolution, and binaries, will be necessary to resolve the issue. 

6. Summary and Conclusions 

Although we have omitted the effects of binaries (which have been shown to retard the 
onset of core collapse, see, e.g., Gao et al. 1991) and also stellar evolution (which can have 
important implications for the early evolution, see Joshi et al. 2000b), the following observations 
and conclusions seem justified: 

• For some two-component clusters the core temperature ratio becomes constant over some 
fraction of the evolution at a minimum value greater than 1, in agreement with previous 
results obtained using other methods. 

• The departure from equipartition calculated for a range of individual and total mass ratios 
approximately respects the theoretical predictions of Spitzer (1969) and Lightman & Fall 
(1978). The agreement with Spitzer is reasonable for m^/m-i > 2, and the Lightman- Fall 
stability boundary (r = 1) appears to reflect the shape of regions of constant £ m i n in 
the parameter space determined by M2/M1 and 1712 /mi, although it lies well inside the 
region occupied by clusters for which equipartition is clearly not attained. A more accurate 
boundary that is suggested by our results is given by equation (|l^). The trend in the 
minimum values predicted for the temperature ratio by Lightman & Fall (1978) is similar to 
what we observe. 

• Stars in the heavier component do not immediately fall into the cluster core as assumed 
by Spitzer in his analysis, and instead many remain well outside the core throughout the 
evolution. 

• The approach to equipartition within the half-mass radius appears to occur on the time scale 
2i cq , as suggested by Spitzer (1969; see eq. (L|). 

• A core temperature ratio of m2/mi appears to be a robust quantity for equilibrium models 
with a relatively shallow potential, in agreement with previous results obtained using other 
methods. 

• For clusters with M2/M1 <C 1 and m^jmx > 1, core collapse times decrease with increasing 
M2/M1 and 7712/771-1, in agreement with previous results. 
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Table 1. Models and Results 



O 


7712/7711 


1VI2 /Mi 


AT 
iV2 


t 


tp/trh 


+ /+ 

<*cc / L rh 


Model 


0.05 


1.50 


0.0272 


1,782 


1.010 


12.8 


12.8 


King 


0.05 


1.75 


0.0216 


1,219 


1.020 


12.0 


12.0 


King 


0.05 


2.00 


0.0176 


876 


1.028 


11.2 


11.2 


King 


0.05 


3.00 


0.00962 


320 


1.027 


9.2 


9.2 


King 


0.05 


4.00 


0.00625 


156 


1.016 


8.7 


9.0 


King 


0.05 


6.00 


0.00340 


57 


1.018 


8.2 


8.2 


King 


0.10 


1.50 


0.0544 


3,502 


1.023 


12.3 


12.3 


King 


0.10 


1.75 


.0432 


2409 


1.035 


10.9 


11.1 


King 


0.10 


2.00 


0.0354 


1,737 


1.043 


9.2 


9.6 


King 


0.10 


3.00 


0.0192 


637 


1.039 


5.9 


6.4 


King 


0.10 


4.00 


0.0125 


311 


1.053 


4.9 


5.5 


King 


0.10 


6.00 


0.00680 


113 


1.061 


3.7 


4.3 


King 


0.15 


1.50 


0.0816 


5,162 


1.027 


11.5 


11.7 


King 


0.15 


1.75 


0.0648 


3570 


1.043 


9.7 


10.1 


King 


0.15 


2.00 


0.0530 


2,583 


1.044 


7.5 


8.5 


King 


0.15 


3.00 


0.0289 


953 


1.082 


4.4 


5.3 


King 


0.15 


4.00 


0.0188 


467 


1.094 


3.3 


3.9 


King 


0.15 


6.00 


0.0102 


170 


1.128 


2.4 


2.6 


King 


0.20 


1.50 


0.109 


6,767 


1.021 


11.0 


11.5 


King 


0.20 


1.75 


0.0864 


4,704 


1.049 


8.4 


9.5 


King 


0.20 


2.00 


0.0707 


3,415 


1.052 


6.8 


7.8 


King 


0.20 


3.00 


0.0385 


1,267 


1.090 


3.8 


4.6 


King 


0.20 


4.00 


0.0250 


621 


1.124 


2.6 


3.2 


King 


0.20 


6.00 


0.0136 


226 


1.138 


1.6 


1.9 


King 


0.25 


1.50 


0.136 


8,318 


1.035 


10.3 


11.0 


Kine: 


0.25 


1.75 


0.108 


5,812 


1.049 


8.1 


9.2 


King 


0.25 


2.00 


0.0883 


4232 


1.072 


6.2 


7.6 


King 


0.25 


3.00 


0.0481 


1,578 


1.104 


3.2 


4.2 


King; 


0.25 


4.00 


0.0313 


775 


1.154 


2.0 


2.6 


King 


0.25 


6.00 


0.0170 


283 


1.175 


1.2 


1.5 


King 


0.50 


1.50 


0.272 


15,358 


1.043 


6.9 


10.2 


King 


0.50 


1.75 


0.216 


10,986 


1.057 


5.0 


8.4 


King 


0.50 


2.00 


0.176 


8,121 


1.079 


3.8 


6.6 


King 


0.50 


3.00 


0.0962 


3,108 


1.173 


2.0 


3.2 


King 


0.50 


4.00 


0.0625 


1,539 


1.312 


1.4 


2.0 


King 


0.50 


6.00 


0.0340 


564 


1.355 


1.0 


1.2 


King 


1.10 


1.50 


0.599 


28,529 


1.046 


1.64 


10.2 


King 


1.24 


5.0 


0.111 


2,172 


1.60 


0.7 


1.2 


King 


0.0042 


2.00 


0.0015 


74 


1.008 


14.5 


14.5 


Plummer 


0.0057 


2.00 


0.0020 


101 


1.021 


14.3 


14.3 


Plummer 


0.0085 


2.00 


0.0030 


150 


1.021 


14.0 


14.0 


Plummer 


0.0101 


2.00 


0.0036 


178 


1.008 


13.9 


13.9 


Plummer 


0.0141 


2.00 


0.0050 


249 


1.021 


13.7 


13.7 


Plummer 


0.0286 


2.00 


0.0101 


503 


1.005 


12.8 


12.8 


Plummer 


0.113 


2.00 


0.040 


1,960 


1.049 


9.0 


9.3 


Plummer 



-18- 



Table 1 — Continued 



s 


milm\ 


M 2 /Mi 


N 2 






^cc/^rh 


Model 


0.314 


2.00 


0.111 


5,262 


1.066 


5.6 


7.5 


Plummer 


2.824 


2.00 


1.000 


33,333 


1.127 


0.014 


7.5 


Plummer 
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Fig. 1. — Evolution of the temperature ratio in the core for S = 0.05 and mijmx = 1.5 (top), 
S = 0.5 and mj/rai = 3 (middle), and S = 1.24 and m-z/mi = 5.0 (bottom). The minimum 
temperature ratio attained in each calculation increases with S. Time is displayed in units of the 
initial half-mass relaxation time (t r h). In each case the evolution is shown until shortly before core 
collapse. Equipartition is clearly not reached prior to core collapse for large S. Notice also that core 
collapse occurs sooner with increasing S. The initial condition in each case was a two-component 
King model with Wo = 6. 
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Fig. 2. — Evolution of the core temperature for the lighter stars (Ti(r < r c ), top) and the heavier 
stars (T-2,(r < r c ), bottom) for the case S = 1.24 and m2/m\ = 5.0. (The ratio of these is shown in 
the bottom-most plot of Fig. |].) Time is displayed in units of the initial half-mass relaxation time 

(trh)- 
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Fig. 3. — Evolution of the core temperature ratio for three calculations with different initial values 
of the dimensionless King model potential Wo, but all with 5 = 1 and m^/m-i = 5 (M2/M1 ~ 0.09). 
From left to right: Wq = 10, Wq = 5, and Wq = 1. While the relaxation and core collapse times for 
these calculations span a wide range, in each case the temperature ratio reaches the same minimum 
value of approximately 1.55. The evolution is shown until shortly before core collapse in each case. 
Note that the logarithmic time scale has compressed the shapes of these curves, so that the leveling 
in the temperature ratio prior to core collapse is not as clearly apparent as in Fig. 1. 
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Fig. 4. — Evolution of the mass density ratio in the core (p C 2/Pci) for the same three cases as in 
Fig. 1: S = 0.05 and m^lmx = 1.5 (top), S = 0.5 and mj/mi = 3 (middle), and S = 1.24 and 
m-ilmx = 5.0 (bottom). The evolution is shown until shortly before core collapse in two cases 
(middle, bottom), and in one case until core collapse (top). As S increases, the time at which equal 
mass densities are attained in the core (t p ) occurs sooner with respect to core collapse (t p ~ t cc for 
the case S = 0.05). 
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Fig. 5. — Evolution of the core mass density for the lighter stars (top) and the heavier stars 
(bottom) for the same case as in Fig. 2: S = 1.24 and m2/m\ = 5.0. (The ratio of these is shown 
in the bottom-most plot of Fig. ||.) The evolution of p c \ is shown until no lighter stars remain in 
the core, whereas the evolution of p& is shown until core collapse. 
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Fig. 6. — Minimum temperature ratio in the core for the 37 calculations in set A, represented here 
using three symbols at points in the parameter space determined by M2/M1 and mil 'm\. Also 
drawn are the Spitzer and Lightman-Fall stability boundaries (S = 0.16 and T = 1, respectively), 
and the boundary A = 0.32 suggested by these results (eq. fl2||). Calculations marked with "□" 
are determined to have reached equipartition within our numerical accuracy. 
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Fig. 7. — Minimum temperature ratio in the core versus S for several values of mil 'mi. Also drawn 
is the Spitzer stability boundary (S = 0.16). 
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Fig. 8. — Minimum temperature ratio in the core versus T for several values of 7712 /mi. Also drawn 
is the Lightman-Fall stability boundary (T = 1) and a theoretical estimate of the minimum core 
temperature ratio, T 1 / 3 ( Lightman fc Fall 1978| ). 
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Fig. 9. — Minimum temperature difference in the core versus the total mass ratio as found in the 
present study (*) and by Inagaki & Wiyanto (1984) (o). Each of these calculations was begun with 
a two-component Plummer model (set B) where mj/mi = 2. Here, vq is the Plummer scale length 
and M is the total cluster mass. Agreement is very good except in the vicinity of small M2/M1. In 
this domain, the temperature of the heavier component is calculated using very few stars, so that 
the difference is characterized by a large amount of noise. 



-28- 



10 



o 



10 



10 



-2 



10 



o 





\ 










+ 







+ 




+ 





























i i i r 



0.10 < tj too < 0.70 



+ 0.70 < tj too < 0.95 



00.95 < t 1 too < 1-00 











+ 
+ 






+ 




+ 





















m 2 /m 1 



+ r=i 



■ 

+ 
o 

_L 



S=0.16 
A=0.32^ 

J I L 



10 



Fig. 10. — Fraction of the core collapse time when equal mass density is attained in the core (t p /t cc ) 
for the 37 calculations in set A, represented here using three symbols at points in the parameter 
space determined by M^/M\ and mj/mi. Also drawn are the Spitzer and Lightman-Fall stability 
boundaries (S = 0.16 and T = 1, respectively), and the boundary suggested by the results shown in 
Fig. |6| (A = 0.32). Note that our boundary appears to apply approximately throughout the range 
1 < m-i/mx < 7 (in contrast to Fig. 6, where it does not apply below mijm\ ~ 1.75). 
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Fig. 11. — Core collapse times versus S for several values of m^/m-i, for 30 calculations in set A. 
The initial condition in each case was a two-component King model with Wq = 6. The times are 
displayed in units of the initial half- mass relaxation time (i r h). 
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Fig. 12. — Minimum core temperature ratio for turning points in a linear series of equilibrium 
models flKatz &: Taff 1983 ). Also drawn is the Spitzer stability boundary (S = 0.16). 
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Fig. 13. — Eight Lagrange radii for the heavier component in the two-component King model with 
S = 1.24, 7712/771-1 = 5 (same model as in Fig. 2). From top to bottom: the radii containing 90%, 
75%, 50%, 25%, 10%, 5%, 1%, and 0.1% of the total mass in the heavier component. Also drawn 
are several points in the evolution of the cluster core radius (•). Note that many stars in the heavier 
component remain well outside of the core throughout the pre-collapse evolution. 



